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Abstract. 

We study certain stationary and time-evolution problems of trapped Bose-Einstein 
condensates using the numerical solution of the Gross-Pitaevskii equation with both 
spherical and axial symmetries. We consider time-evolution problems initiated by 
changing the interatomic scattering length or harmonic trapping potential suddenly in 
a stationary condensate. These changes introduce oscillations in the condensate which 
are studied in detail. We use a time iterative split-step method for the solution of the 
time-dependent Gross-Pitaevskii equation, where all nonlinear and linear nonderivative 
terms are treated separately from the time propagation with the kinetic energy terms. 
Even for an arbitrarily strong nonlinear term this leads to extremely accurate and 
stable results after millions of time iterations of the original equation. 
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1. Introduction 

Since the successful detection [|1| of Bose-Einstein condensates in dilute weakly- 
interacting bosonic atoms employing magnetic trap at ultra-low temperature, there 
have been intense theoretical studies on different aspects of the condensate [^-|i"3| . The 



experimental magnetic trap could be either spherically symmetric or axially symmetric. 
The properties of an ideal condensate at zero temperature are usually described by 
the time-dependent, nonlinear, mean-field Gross-Pitaevskii (GP) equation jl4| which 



incorporates appropriately the trap potential as well as the interaction between the 
atoms forming the condensate. There are many numerical methods for the solution of 
the GP equation fl-lOU- The effect of the interatomic interaction leads to the nonlinear 



term in the GP equation which complicates the solution procedure. Also, to simulate the 
proper experimental situation one should be prepared to deal with an axially-symmetric 
harmonic oscillator trap in addition to a spherically symmetric one [ff|-[T0|]. 

A numerical study of the time-dependent GP equation is extremely interesting from 
a physical point of view, as this can provide solution to many stationary as well as time- 
evolution problems involving the condensate. The stationary problems are governed by 
a wave function with trivial time (r) dependence \&(r, r) = ^(r) exp(— ifir/h), so that 
|\]/(r, r)| = |\&(r)|, where \i is a real energy parameter. Such problems with trivial time 
dependence and under the action of a spherically symmetric trap can be treated by 
time-independent Runge-Kutta integration method ||. The solution of these problems 
can also be extracted from the solution of the time-dependent GP equation In the 

presence of an axially-symmetric trap, there is a time-independent expansion scheme 
for the solution of the stationary problem ||. However, for most applications one solves 
the time-dependent GP equation and extract the stationary wave function \l/(r) and 
the parametric energy // [p|,[10|. In addition to obtaining the solution of the stationary 
problem, the time- dependent GP equation can be used to study the intrinsic time- 
evolution problems with nontrivial time dependence. 

There are several methods for the numerical solution of the time-dependent GP 
equation. Most of them adopt a similar time-iteration procedure in execution, although 
the details could be different. The time-dependent GP equation is first discretized in 
space and time and then solved iteratively with an initial input solution j|,^,|9],[10| . In the 



spherically symmetric case the GP equation is solved by discretization with the Crank- 



Nicholson scheme [|I7^|T6|j for time propagation complimented by the known boundary 
conditions ||. However, in the axially symmetric case a two-step procedure |TjJ is 
first used to separate the radial and axial parts of the Hamiltonian before applying 
the the explicit J7J or the Crank-Nicholson scheme [|3|-|T0[| for time propagation in each 
step. In all these approaches the nonlinear term is treated together with the explicit or 
the Crank-Nicholson scheme P, jo], |5|-|TT)|j . It is well known that such a Crank-Nicholson 
method may suffer from a possible numerical amplification of random numerical noise, 
due to nonlinearity, leading to spurious input of kinetic energy and limiting the stability 
of the algorithm to short time (small number of iterations) |§. This also limits the 



Bose-Einstein condensation dynamics 



3 



applicability of the method only to small values of the nonlinearity ||]. 

In this paper we suggest a different split-step time-iteration method for the solution 
of the time-dependent GP equation where the full Hamiltonian is split into three 
parts [injj . The nonlinear as well as the different linear (nonderivative) terms (excluding 
spatial derivatives) are treated in one step. In the spherical case the spatial derivative is 
treated in another step. In the axially symmetric case the radial and axial derivatives are 
dealt with in two separate steps. In the treatment of nonderivative terms the numerical 
error is kept to a minimum. In dealing with the derivative kinetic energy terms the 
GP equation is solved by discretization using the Crank-Nicholson scheme for time 
propagation complimented by the known boundary conditions [ 10|, 16j ,|T7 1 . The advantage 
of the present split-step procedure is that the nonlinear and other linear nonderivative 
terms can be treated very precisely and this improves the accuracy and the stability of 
the method compared to previous ones [[UJ. Consequently, we could easily find accurate 
solution of the GP equation for very large nonlinearity, which remains stable even after 
millions of iterations. This will be of advantage in using the present approach for the 
study of time-evolution problems during a large interval of time. 

As an application of the present approach in studying time-evolution problems we 
consider in the spherically symmetric case the following two types of problems: (1) On 
a previously formed condensate the harmonic oscillator trapping potential is increased 
or decreased suddenly by a factor of two. This can be achieved in the laboratory 
by changing the current in the coil responsible for the magnetic field responsible for 
trapping. (2) On a previously formed condensate the atomic scattering length is 
increased or decreased suddenly by a factor of two. Case (2) above is of interest as 
recently it has been possible to modify the scattering length in a controlled fashion by 
exploiting a Feshbach resonance |T8| , [I5|] . In both cases we study the resultant oscillation 
of the root mean square radius of the condensate and its first time derivative. In the 
axially symmetric case we also consider the time-evolution problems above and study 
the resultant oscillation of the root mean square sizes in radial and axial directions. 

In section 2 we describe briefly the time-dependent GP equation with spherical and 
axial traps. The split-step Crank-Nicholson method for time propagation is described in 
section 3. In sections 4 and 5 we report the numerical results of the present investigation 
for the spherically and axially symmetric cases, respectively. Finally, in sections 6 and 
7 we present a discussion and summary of our study. 



2. Nonlinear Gross-Pitaevskii Equation 



At zero temperature, the time- dependent Bose-Einstein condensate wave function 
^(r; r) at position r and time r may be described by the following mean- field nonlinear 
GP equation P,|H 



2m 



+ V{r)+gN\^(r;T) 



0. 



(2.1) 
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Here m is the mass and iV the number of atoms in the condensate, g = AnTi 2 a/m 
the strength of interatomic interaction, with a the atomic scattering length. The 
normalization condition of the wave function is J dr|\I/(r; r)| 2 = 1. 



2.1. Spherically Symmetric Case 

In this case the trap potential is given by V^(r) = ^mu 2 r 2 , where u is the angular 
frequency and r the radial distance. The wave function can be written as ^/(r;r) = 

i/j(r,r). After a transformation of variables to dimensionless quantities defined by 

x = V2r/l, t = tu, I = J {Ti/mui) and (p(x,t) = (p(x,t)/x = ip(r, r)(47r/ 3 ) 1 / 2 , the 
GP equation in this case becomes 

„2 



o- 



dx 2 



x 

T 



ip(x,t) 


2 d~ 


X 





(p(x,t) = 0, 



where M = Na/l. The normalization condition for the wave function is 
dx\(f{x,t)\ 2 = 2y/2. 



(2.2) 



(2.3) 



2.2. Axially Symmetric Case 

The trap potential is given by V(r) = \muj 2 (p 2 + A 2 z 2 ) where u is the angular frequency 
in the radial direction p and Xuj that in the axial direction z. We are using the 
cylindrical coordinate system r = (p, 6, z) with 9 the azimuthal angle. In this case 
one can have quantized vortex states with rotational motion around the z axis. In 
such a vortex the atoms flow with tangential velocity TiL/(mr) such that each atom has 
quantized angular momentum UL along z axis. The wave function can then be written 
as: \&(r; r) = ip(r, z; r) exp(iL60, with L = 0, ±1, ±2, .... 

Using the above angular distribution of the wave function in (|2.1|), in terms 
of dimensionless variables x = \/2r/l, y = \/2z/l, t = tu, I = yh/ (mu), and 
cj)(x,y;t) = (f(x,y;t)/x = y/4irl 3 tp(r, z; r), we get 

d 2 Id d 2 



ld_ 

dx 2 x dx 



dy 2 



~ + l{ X +XlJ ~~ 



X 



N 



ip{x,y;t) 



x 



.d_ 

'di 



<p(x,y,t) = 0. 



The normalization condition of the wave function is 

dx \ dy\tp{x,y\£)fx- x = 2V2. 
Jo 



(2.4) 



(2.5) 



3. Split-Step Crank-Nicholson Method 



The GP equation in both spherically symmetric and axially symmetric cases can be 
formally written as 



■9<P TT 



(3.1) 
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where the Hamiltonian H contains the different nonlinear and linear terms including 
the spatial derivatives. We solve this equation by iteration [H^-jlT]. A given trial input 
solution is propagated in time over small time steps until a stable final solution is 
reached. The GP equation is discretized in space and time using the finite difference 
scheme. This procedure results in a set of algebraic equation which can be solved by 
time iteration using an input solution consistent with the known boundary condition. 
In the present split-step method this iteration is conveniently done in several steps 
by breaking up the full Hamiltonian into different derivative and non-derivative parts. 

In the spherically symmetric case we split H in to three parts: H = Hi + H 2 + H 3 , 
where 



2 





<p(x,t) 


2" 




X 





H 2 
H 3 



d 2 



dx 2 



Hi. 



(3.2) 

(3.3) 
(3.4) 



The time variable is discretized as t n = nA where A is the time step. The solution is 
advanced first over the time step A at time t n by solving the GP equation ( |3.1| ) with 
H = Hi to produce an intermediate solution (p n+1 / 3 from Lp n , where tp n is the discretized 
wave function at time t n . As there is no derivative in Hi this propagation is performed 
essentially exactly for small A through the operation 

^n+i/3 = O nd (Hi)tp n = e- iA *V, (3.5) 

where O n ^(Hx) denotes time-evolution operation with Hi and the suffix 'nd' denotes 
non- derivative. Next we perform the time propagation corresponding to the operator 



H 2 numerically by the following semi-implicit Crank-Nicholson scheme |L6 

n+2/3 _ ^n+1/3 j 

~ ~ 2 



-iA 



-H 2 (ip n+2 ^ + ip n+1 / 3 ). 



The formal solution to ( p.6| ) is 

^ n+2/3 = <?CN(^V n+1/3 



1 - iAH 2 /2 1/8 

) r ) 



(3.6) 



(3.7) 



1 + iAH 2 /2 

where Oqj^ denotes time-evolution operation with H 2 and the suffix 'CN' refers to 
the Crank-Nicholson algorithm. Operation Oq-^ is used to propagate the intermediate 
solution (p n+1 / 3 by time step A to generate the second intermediate solution Lp n+2 / 3 . The 
final solution is obtained from 

V n+l = O^mO^mO^Hi)^. (3.8) 

The break-up of the nonderivative term in two parts —Hi and -£^3— symmetrically 
around the derivative term H 2 , increases enormously the stability of the method and 
reduces the numerical error. 
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In the axially symmetric case H is split into three parts: H = Hi + H 2 + H 3 , where 

rr 1 ( 2 . X 2 2 4 \ Z, 2 <p(x,y,t) 2 , , 

H ' = A X +Xy + —r~ ' (3 - 9) 

d 2 1 9 

H 2 =-^ + -—, (3.10) 

(3 ' n) 

The strategy is now obvious. The solution is advanced first over the time step A at 
time t n using the Hamiltonian H\ above to produce an intermediate solution y? n+1 / 3 
from ip n . For small A this propagation is performed essentially exactly via ( |3.5|) . Next 
we perform two successive time propagations corresponding to the operators H 2 and H 3 
numerically by the Crank-Nicholson scheme (|3.7|) . Consequently, a single time iteration 
from t n to t n+ \ is performed via the following three-step operation: 

^ n+1 = 0CN(#3)0CN(^) nd(tfiV n - ( 3 - 12 ) 
In the axially symmetric case due to the intrinsic asymmetry of the kinetic energy terms 
in x and y directions, the complete symmetrization as in the spherically symmetric case 
( P-8| ) is nontrivial and we did not attempt such a symmetrization. 

The advantage of the above split-step method with small time step A is due to 
the following three factors [15,|l6|. First, all iterations conserve normalization of the 
wave function. Second, the error involved in splitting the Hamiltonian is proportional 
to A 2 and can be neglected and the method preserves the symplectic structure of the 
Hamiltonian formulation. Finally, as a major part of the Hamiltonian including the 
nonlinear term is treated fairly accurately without mixing with the Crank-Nicholson 
propagation, the method can deal with an arbitrarily large nonlinear term and lead to 
stable and accurate converged result. 

Now we describe the Crank-Nicholson algorithm in the spherically- and axially- 
symmetric cases. In the spherically symmetric case the GP equation is mapped on to 
N x one-dimensional grid points in x. Equation (|3.1| ) is discretized with H = H 2 of ( |3.3j ) 
by the following Crank-Nicholson scheme [JI31-0]: 

^^T 1 - <$) _ 



A 2h 2 



/,„n+l o,„n+l i ,„n+l\ i / , _n o,^" i ,„W 



(3.13) 



where <£>™ = <p(xj,t n ) refers to x — Xj = jh, j = 1, 2, N x and h is the space step. This 
scheme is constructed by approximating d/dt by a two-point formula and d 2 /dx 2 by a 
three-point formula. This procedure results in a series of tridiagonal sets of equations 
( p. 13 ) in ipj+l, <fi] +1 , and ipjll at time t n+ i, which are solved using the proper boundary 
conditions. 

The Crank-Nicholson scheme ( |3.13|) possesses certain properties worth mentioning 



T^,|16l. This error in this scheme is both second order in space and time steps so that 



for small A and h the error is negligible. This scheme is also unconditionally stable. 
The boundary condition at infinity is preserved for small values of A/h 2 [ |16|1 . 
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In the axially symmetric case the time-dependent GP equation ( |2.4j ) is mapped on 
a two-dimensional grid of points N x x N y in x and y. Equation (|3.1| ) with H = H 2 oi 
( p.lOj ) is discretized using the following Crank-Nicholson scheme [15-17]: 

A 2h 2 



+ 



Axjh 



(3.14) 



where the discretized wave function tp^ p = (p(xj,y p ,t n ) refers to a fixed y = y p = ph, 



P = l,2, 



iV y at different x 



j/i, j = 1,2, ...,N X , and is the space step. The 



error in scheme (|3.14) is also second order in both space and time steps. This procedure 



results in a tridiagonal set of equations ( |3.14| ) in y?"+i iP , , and Pjli p at time t n+ \ 



for each y p , which are solved consistent with the boundary conditions [17]. Equation 
Q3.1| ) with H = H3 of ( 3 . 1 1|) is discretized similarly: 

? Y(o" +1 - m n ) 1 r 1 

A = ~2W ~ Vj ' p + V ^ + ~ + V ^-i>\ ' ( 3 - 15 ) 

where now ipj refers to a fixed Xj = jh for all y p = ph. Using the solution obtained after 
x iteration as input, the discretized tridiagonal equations ( j3.15|) along the y direction 
for constant x are solved similarly. 

The iteration is started with the following normalized analytic solution corres- 
ponding to the ground states of ([2.2|) and (|2.4|) for spherical and axial symmetry, 



respectively, with the coefficient of the nonlinear term set to zero 



tp(x) = n-°- 25 2xexp(-x 2 /A) 



and 



A 



2 2 I L I(|L|!) 2 



7T 



0-25 x l+| L | 



27T 



exp[-(x 2 + Ay 2 )/4]. 



(3.16) 



(3.17) 



The GP equation was discretized with space step h — 0.1 and time step A = 0.001. For 
small nonlinearity the largest values of x and y are x ma x = 8, |y|max = 8- However, 
for stronger nonlinearity, larger values of x max and |y|max (up to 25) are employed. 
The norm of the wave function is conserved after each iteration due to the unitarity of 
the time evolution operator. However, it is of advantage to reinforce numerically the 
proper normalization of the wave function given by (|2.3|) or ( |2.5|) after each complete 
time iteration in order to improve the precision of the result. During the iteration the 
coefficient of the nonlinear term is increased from at each step by A x = 0.0001 for 
both spherical and axial symmetry until the final value of nonlinearity M is attained at 
a time called time t = 0. This corresponds to the final solution. Then about a million 
time iterations of the equation were performed which shows the stability of the result 
at t — 1000. The numerical values of the different steps and parameters were fixed after 
some experimentation in order to reduce the numerical error. 

The maximum fluctuation of the result was found to be below few percent in the 
course of a million iterations. The final wave function was found to be very smooth 
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at a specific time in the course of time iteration. Consistent with the dynamics, small 
and steady oscillation of the result appears as the nonlinear term is increased by step 
Ai = 0.0001. This dynamical oscillation mostly accounts for the error in the result, 
which can be reduced by reducing the parameter A x . We verified that there was no 
increase of this numerical error with time iteration. 

For large nonlinearity, the Thomas-Fermi (TF) solution of the GP equation is a 
better approximation to the exact result || than the harmonic oscillator solutions (|3.16| ) 
and ( gl7p . In that case it is advisable to use the TF solution as the initial trial input 
to the GP equation with full nonlinearity and consider time iteration of this equation 
without changing the nonlinearity |J. This time iteration is to be continued until a 
converged solution is obtained. However, in all the calculations reported in this paper 
only Q3.16Q and (|3.17|) are used as trial inputs. 



4. Result for Spherically Symmetric Trap 

4-1- Stationary Problem 

First, we study the stationary problem in the spherically symmetric case. In figure 1 (a) 
we plot the wave function (p(x, t) for four different values of nonlinearity M = 5, 100, 1000 
and 5000. The plotted wave functions are normalized according to / °° dxx 2 \(j)(x, t)\ 2 = 
2y / 2A/' and not according to (|2.3|) . The corresponding statistical error in the wave 
function is plotted in figure 1 (b). The error was calculated through the course of 10 6 
time iterations with time step A = 0.001 of the final solution generating a time interval 
of 1000 units of time t. A sample of the wave function was chosen after 1000 iterations 
each. The percentage statistical error E was calculated from these samples of wave 
function via 



E 



100 

<f)(x) 



j 



7 E (l<M*)| 2 -0(*) 

J i=i 



1/2 



(4.1) 



where <f) is the average of \<p\ over total number of samples J = 1000. From figure 1 
(b) we find that the percentage error is less than about 0.5% for the central part of 
the condensate where the wave function is sizable. It increases to about 2% near the 
boundary where the wave function is almost zero. The most interesting feature of the 
present calculation is that the solution remains stable even after a million iterations. 
We calculated the statistical error for the first and last 100 of the 1000 samples above. 
The errors so calculated were the same as those in figure 1 (b), which shows that error 
does not increase with time. In other previous approaches it was not possible to obtain 
stable result after such a large number of iterations of the GP equations |10| , [TT[| . The 
stability with time iteration in the present approach will be an asset in the study of 
time evolution problems which we undertake next. 



Bose-Einstein condensation dynamics 



9 



T 1 1 1 1 1 1 1 1 1 1 1 1 1 1 - 
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n=5000 




j i i i i i i i i i i i i i i_ 



Figure 1. The wave function in the spherically symmetric case (a) \4>(x)\ vs. x at 
time t — 1000 for nonlinearity AT = 5, 100, 1000, 5000 after a million time iterations of 
step A = 0.001 of the solution at t — and (b) the corresponding average statistical 
error vs. x during this interval. The plots are labeled by the respective n values and 
the wave functions are normalized as / °° dxx 2 \4>(x)\ 2 — 2\f2M . 



4-2. Time Evolution Problems 

There are several interesting time evolution studies that one can pursue. After the 
generation of a stable condensate one can reduce or increase suddenly the strength of 
the harmonic oscillator trap by a factor and study the oscillation of the condensate 
thereafter. Also, now it is possible to manipulate the effective strength of interatomic 
interaction through a Feshbach resonance after a stable condensate has been formed by 
changing the surrounding electromagnetic field 18|, 19|. So one can increase or reduce 
suddenly the strength of the nonlinear term through a change of the scattering length 
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in this fashion and study the oscillation of the condensate thereafter. In both these 
cases one observes the time evolution of the root mean square (rms) radius X of the 
condensate. This steady oscillation is best studied theoretically by plotting 'velocity' X 
vs. 'position' X, where dot denotes time derivative in the course of time evolution. 

3 
2 
1 

■X 
-1 
-2 
-3 

6 7 8 9 10 

X 

Figure 2. Root mean square velocity X vs. X during sustained regular periodic 
oscillation of the spherically symmetric condensate with AT = 2000 initiated at t = 
by suddenly (a) increasing or (b) decreasing the harmonic oscillator term by a factor 
of 2, and by suddenly (c) increasing or (b) decreasing the nonlinear term by a factor 
of 2. 

For the time evolution study we consider a previously formed condensate with 
M = 2000 prepared at t = as in the study of the stationary problem. We then inflict 
the four following changes in the system. At t = we (a) increase or (b) decrease 
suddenly the coefficient of the harmonic oscillator x 2 /4 term in (|2.2|) by a factor of 2. 
Next at t = we (c) increase or (d) decrease suddenly the coefficient of the nonlinear 
term M in (|2.2j ) by a factor of 2. In all these four cases we iterate the GP equation in 
time with time step A = 0.001 and observe the system for an interval of 1000 units of 
time over 10 6 iterations. The maximum radial distance was taken to be £max — 25. 

In figure 2 we plot velocity X vs. radius X in the above four cases. We obtain 
a closed curve in the phase space which confirms the periodic nature of a sustained 
oscillation of the condensate over 1000 units of time. The clean closed curve over such 
a long interval of time assures of the low numerical error in the calculation. When the 
harmonic oscillator term is doubled or the nonlinearity halved, the system is compressed 
and the rms radius oscillates between its initial value and a smaller final value. When 
the harmonic oscillator term is halved or the nonlinearity doubled, the system expands 
and the rms radius oscillates between its initial value and a larger final value. The four 
curves of figure 2 should originate at a central point as the initial condition on all of 
them is at t = 0, X = and X = constant. However, the intersection of the pair 
of curves at the other extremities is accidental. This periodic oscillation of the wave 
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Figure 3. The wave function |0(a;,t)| of the spherically symmetric condensate vs. x 
and t for AT = 2000 during breathing oscillation in case of figure 2 (a). 

function is clear from the phase space diagram presented in figure 2. This oscillation 
is explicitly shown in figure 3 where we plot \<f>(x, t)\ vs. x and t in case (a) above 
when the harmonic oscillator term is doubled suddenly at t = on a previously formed 
condensate with M = 2000. The regular breathing-mode-type pattern of oscillation of 
figure 3 apparently continues for ever. The periodic increase and decrease in central 
density are accompanied, respectively, by a periodic decrease and increase in radius. 
Similar oscillation exists in the three other cases considered in figure 2 above. 

We also calculated the frequency of oscillation of the rms radius X in the four cases 
considered in figure 2. The frequencies of simulation in cases (a), (b), (c), and (d) above 
are 0.50, 0.25, 0.35, and 0.35. We are measuring time in units of oj~ x or (2-ku)^ 1 , where 
v is the frequency of the harmonic trapping potential. Hence in present units with u = 1 
the trap frequency corresponding to (|2.2|) is v — \/{2it). We studied similar oscillation 
as in cases (a) and (b) above for nonlinearity M = 0. In both cases the frequency of 
oscillation was two times the existing harmonic oscillator frequency. As the frequencies 
are increased and reduced by y/2 in cases (a) and (b), two times the existing frequency 
for (a) is 2u = \f2jix ~ 0.45 and for (b) is 2v = l/(V2n) ~ 0.23. In case of (c) and (d) 
the frequency is unchanged and 2v = 1/tt ~ 0.32. These numbers compare well with 
the respective results of simulation, e.g., 0.50, 0.25, and 0.35. The difference between 
the two sets is due to the large nonlinearity (jV = 2000) present in simulation. 

5. Result for Axially Symmetric Trap 

5.1. Stationary Problem 

We study the stationary problem in the axially symmetric case. All results reported in 
this section are calculated with A = a/8. For L = 0, we calculated the wave function 
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4>(x,y) for nonlinearities up to TV = 500 and normalized according to ( |2.5| ). We also 
calculated the percentage statistical error over 100 units of time calculated from 100 
samples via ( |4.1|) . For all values of nonlinearities, the error is found to be smaller than 
few percentage points almost over the entire region. As in the spherically symmetric 
case, the error in the central region remained much less than near the peripheries. Even 
for the largest nonlinearity the convergence is smooth and the error remained small. 




Figure 4. The vortex wave function in the axially symmetric case (a) \<fi(x,y)\ vs. 
x and y at time t — 100 for M = 100 and L = 2 after 100000 time iterations of 
step A = 0.001 after forming the solution at t = and (b) the corresponding average 
statistical error vs. x and y during this interval. The wave function is normalized 
according to (2.5). 



We also tested our approach for vortex states, which could be more difficult to 
calculate numerically due to the 1/x 2 term. In figure 4 (a) we plot the vortex wave 
function for L = 2 and Af = 100. The corresponding numerical error calculated with 
100 samples over 100 units of time is plotted in figure 4 (b). Even for the vortex state 
L = 2, with a nonlinearity as large as M = 100 the error is very small as we can see in 
figure 4. 



5.2. Time Evolution Problem 

Next we study time evolution problems with an axially symmetric trap. We take an 
initial state with J\f = 50 and A = which is prepared as usual at time t = 
when we reduce suddenly the strength of the nonlinear term by a factor of 2. This 
corresponds to reducing the scattering length by the same factor, which can be realized 
in experiments []18|)[19|. The system then starts to oscillate which may continue for ever 
as in the spherically symmetric case, albeit with different frequencies in radial and axial 
directions. We calculate the time evolution of the rms values of x and y: X and Y. The 
oscillation is more involved and in figure 5 we plot X and Y vs. t. Finally, we consider 
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again the same initial state with M = 50 at t = and we double suddenly both the 
radial and axial harmonic oscillator trapping terms. The system starts to oscillate after 
the change. This is clearly exhibited by plotting X and Y vs. t in figure 6. 



2.5Y 




, , , , i , , , , i , , , , i , , , , i , , , , 

300 305 310 315 320 325 

t 



Figure 5. Root mean square sizes X and Y vs. t during the oscillation of a previously 
formed axially symmetric condensate with Af = 50 when the nonlinearity is suddenly 
reduced by a factor of 2. 



We calculated the frequency of oscillation of rms sizes X and Y. In case of the 
simulation corresponding to figure 5 the principal frequency of oscillation along X 
direction is V\ ~ 0.29 and that in Y direction is v% ~ 0.84. In case of simulation of 
figure 6 these frequencies are v\ ~ 0.42 and v 2 ~ 1.16, respectively. We also considered 
the frequency spectrum of the time variations represented in figures 5 and 6 and we noted 
se veral freq uencies, of which the principal in each case are, v\, u 2) 2u\, 2z/ 2 , \u 2 ± v\\ and 

Similar oscillation has been observed in a recent experiment [19| for a dynamically 
changing axially symmetric condensate with small nonlinearity — the observed 
frequencies being twice the existing harmonic oscillator trap frequencies ^radial and 
^axial a l° n S radial and axial directions. We repeated the simulation of figure 6 with 
nonlinearity zero and find that the resultant frequencies of simulation are exactly double 
the existing harmonic oscillator trap frequencies ^ rac ii a i and ^ ax j a i along radial and 
axial directions. In case of figure 5 these frequencies are 2i/ ra( ij a j = 1/n ~ 0.32 and 
^axial = Vo/7r ~ 0.90 (^ ax i a i/ ^radial = ~ 2.83), which are to be compared 
with the frequencies V\ ~ 0.29 and u 2 ~ 0.84 {v^/vi ~ 2.90) obtained in the present 
simulation, respectively. The small disagreement between the two sets is due to the 
large nonlinearity of the present simulation. In case of figure 6 the existing trap 
frequencies are increased by a/2 and hence the expected frequencies of oscillation are 
^^radial = V^/ 71 " ~ 0.45 and 2z/ ax j a ^ = 4/tt w 1.27, which compare well with the 
following results of simulation with large nonlinearity: V\ ~ 0.42 and v 2 ~ 1.16, 
respectively (z^/^i ~ 2.76), again the difference being due to the strong nonlinearity. 
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Figure 6. Root mean square sizes X and Y vs. t during the oscillation of a previously 
formed axially symmetric condensate with TV = 50 when the harmonic oscillator terms 
are suddenly increased by a factor of 2. 

Although, we could not relate the present frequencies of numerical simulation to those 
of the low-energy elementary excitations of a Bose-Einstein condensate in an axially 
symmetric trap ||20|| , it is interesting to note that the nonlinear combination J {y\ + u$) 
appears in both the cases. 




6. Discussion 



In this paper we have presented an account of the split-step Crank-Nicholson method 
for the solution of the GP equation, where the nonlinear and all nonderivative terms are 
treated in a separate step essentially exactly via (|3.5|). In previous applications P,[5|,|7|-4TC|] 
of the Crank-Nicholson method to the GP equation these terms were included in the 
finite-difference scheme for the space derivatives. In the pioneering application of this 
approach Ruprecht et al 0] noted that the presence of the nonlinear potential sets a limit 
to the numerical stability of the method and limits the ability to find the wave function 
for large nonlinearity. For the spherically symmetric case, Ruprecht et al |3| produced 
results for values of nonlinearity J\f = Na/l upto 25 y2. For the axially symmetric 
case, the maximum of M employed by Dalfovo and Modugno || and by Holland and 
Cooper H were 20 and 5, respectively. Using a different method Chiofalo et al [[0| 



produced results for values of M upto 6. However, unlike in these other methods, in the 
present method the results remain stable for large A/", 

At this point we must comment on the potentially powerful split-operator method 
with a pseudo-spectral scheme for derivatives ||21|| , often used in solving nonlinear 
Schrodinger-type equations. In the pseudo-spectral scheme the unknown function is 
expanded accurately in terms of some known polynomial (e. g. Chebyshev polynomial) 
over a set of grid points. As the number of terms in this expansion increases, the discrete 
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representation of the function and its derivatives becomes increasingly accurate. In the 
present Crank-Nicholson method the space derivative is approximated by a simple three- 
point finite-difference scheme. To the best of our knowledge so far there has been no 
systematic application of the pseudo-spectral method to the solution of the GP equation, 
specially for the axially symmetric case |22| . Whether the added numerical effort in the 



pseudo-spectral method is compensated by the increased accuracy can only be found 
after a systematic comparison of the two approaches. This remains a problem of future 
interest . 

7. Summary 

In this paper we propose and implement a split-step method for the numerical solution of 
the time-dependent nonlinear GP equation under the action of a trap with both spherical 
and cylindrical symmetries by time propagation starting with the initial input solution 
for the harmonic oscillator problem. The full Hamiltonian is split into the derivative 
and nonderivative parts. In this fashion the time propagation with the nonderivative 
parts can be treated very accurately. The spatial derivative parts are treated by the 
Crank-Nicholson method. In the axially symmetric case the spatial derivatives in the 
radial and axial directions are dealt with in two independent steps. This split-step 
method leads to highly stable and accurate results. The final result remains stable for 
millions of time iteration of the GP equation. 

We applied the above method for the numerical study of certain stationary and 
time-evolution problems with spherical and axial traps. In the spherical case stationary 
solutions with nonlinearity Af = Na/l = 5, 100, 1000, and 5000 showed accurate result 
with small numerical error. For a condensate with Af = 2000 we studied certain time- 
evolution problems. We studied the resulting oscillation of the condensate when at 
t = the trap frequency or the scattering length is altered suddenly. In all clean 
periodic motion of the root mean square radius was noted. 

In the axially symmetric case the stationary solutions for L = with nonlinearity 
Af = Na/l = 50 and 300 showed accurate result. A calculation of the vortex state 
for L = 2 with nonlinearity Af = 100 showed equally precise and stable results. The 
numerical error was less than one percent in the central region and a couple of percents 
in the peripheral region. We also performed two time-evolution studies for Af = 50 when 
the harmonic oscillator trap is doubled or the scattering length is reduced suddenly. The 
root mean square sizes along radial and axial directions showed clean periodic oscillation 
which is studied in detail. The present method seems to be very attractive for studying 
various dynamical evolution problems with dissipative (imaginary) interaction, which 
can not be handled efficiently by other methods such as variational and Thomas-Fermi 
methods. We have already made such an application in the study of chaotic dynamics 
of the Bose-Einstein condensate using the Gross-Pitaevskii equation with dissipative 
interaction f23| . 
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